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Abstract 

The space-uniform amplitude envelope of the Ion Temperature Gra- 
dient driven turbulence is unstable to small perturbations and evolves to 
nonuniform, soliton-like modulated profiles. The induced poloidal asym- 
metry of the transport fluxes can generate spontaneous poloidal spin-up 
of the tokamak plasma. 
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1 Introduction 

The ion temperature gradient turbulence, considered a major source of anoma- 
lous transport in the tokamak plasma, is characterized by the coexistence of 
irregular patterns (randomly fluctuating field) and intermittent robust cuasi- 
coherent structures. In closely related fluid models (for exemple, in the physics 
of atmosphere) modulational instabilities are known to produce solitary struc- 
tures on the envelope of the fluctuating field. In the case of the tokamak, this 
can be particularly important since a poloidally nonuniform amplitude of the 
turbulence generates nonuniform transport rates. For a sufficiently high nonuni- 
formity the torque arising via the mechanism initially mentioned by Stringer jlj 
may overcame the rotation damping due to the poloidal magnetic pumping. In 
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this work we show that this can indeed be the case by proving that the uni- 
form poloidal envelope of the ITG turbulence is an unstable state. It appears 
that the ITG turbulence has intrinsic resources to generate poloidal rotation 
via a combination of mechanisms which is not directly related to the Reynolds 
stress, inverse cascade, direct-ion loss, or other classical sources of rotation. We 
provide an essentially mathematical explanation of the instability of the solu- 
tion consisting of poloidally uniform envelope of the turbulence. Based on the 
geometrico-algebraic method of solving integrablc nonlinear differential equa- 
tions on periodic domains we invoke an existing result, that any perturbation 
removes the degeneracies due to coincident eigenvalues in the main spectrum of 
the Lax operator, thus changing the topology of the hyperelliptic Riemann sur- 
face that provides the solution. The perturbed solution separates exponentially 
(in function space) from the initial one (uniform envelope) and this yields an 
exponential growth of the poloidal nonuniformity. The magnitude of the torque 
may be comparable to the poloidal magnetic damping. 

There are many works in plasma physics related to the soliton dynamics Q , 
[|| , Q , || , , etc. In a recent paper the focusing solution has been used to study 
the formation of coherent motion or intermittent patterns (streamers) (^]. Our 
work is basically a mathematical approach, using well developed technics related 
to the Inverse Scattering Transform method for periodic domains. However, 
more physical analysis may become possible combining the knowledge of the 
spectral properties of the ion turbulence with the nonlinear stability properties. 

This work consists of several parts that may appear as developing separately: 
the derivation of the ion equation (barotropic equation), the multiple space- 
time scales analysis, the solution of the Nonlinear Schrodinger Equation and 
the stability of the solution, the torque arising from the Stringer mechanism 
and the possibility of the rotation. For being self-contained the paper also 
includes a review of the multiple space-time analysis (from the atmospheric 
physics applications) and of the geometric-algebraic method of integration of 
the Nonlinear Schrodinger Equation. Even if these parts could be found in basic 
works (see the references) they are included both for clarity and for their extreme 
importance for further applications related to our problem or independent of 
this. 

2 The slab model of the ion mode instability 

We consider cylindrical geometry with circular magnetic surfaces. Locally the 
model can be reduced to a slab geometry with (x, y) cartezian coordinates re- 
placing respectively the radius and poloidal angle coordinates (r, 9). At equilib- 
rium the plasma parameters are constant on the magnetic surfaces. The effects 
of the toroidicity and of the particle drifts are not included instead the nonlin- 
earity related to the ion polarization drift is fully retained. The plasma model 
consists of the continuity equations and the equations of motion for electrons 
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and ions: 



^ + V- (n Q v Q ) = (1) 

m a n a + (v Q • V) = -Vp Q - V • n a + e a n a (E + v Q xB) + R Q 

(2) 

where a — e,i. The friction forces R e = — Hi = —n\e\J\\/a\\ , which arc im- 
portant for the parallel electron momentum balance, vanish for infinite plasma 
conductivity, which we will assume. The collisional viscosity ir e j will be ne- 
glected as well. However we will need to include it later when we will consider 
the balance of the forces contributing to the poloidal rotation. The electron and 
ion temperatures are considered constant along the magnetic lines Vi|T ei , = 0. 
The equilibrium quantities are perturbed by the wave potential <f>: n = no + n. 
A sheared poloidal plasma rotation is included, and we later will make explicit 
the corresponding part in the potential, <j> . 

The momentum conservation equations are used to determine the perpendic- 
ular velocities of the electrons and ions. The parallel momentum conservation 
equation for electrons, in the absence of dissipation or drifts gives the adiabatic 
distribution of the density fluctuation. The velocities are introduced in the con- 
tinuity equations to find the dynamical equations for the density and electric 
potential. 

From the equations of motion for the ions the velocities are obtained in the 
form: 

v» = v_u = Vdia,i + v E + v poM 
where the ion diamagnetic velocity is 

1 1 - TT 
Vdia,i = 7-11 X Vp t 

Hi rriiUi 

The versor of the magnetic field is n; the versors along the transversal coordinate 
axis (x,y) will be noted (e x ,e y ). The ion-polarization velocity is: 

v poM = -^w(^ + (v E 'V)v,)xn (3) 
d \ 



BCli \dt 

Using the notation ^ e ^ + (vg • Vj_) the perpendicular ion velocity can be 
written 

-V<j> x n 

V_L,i = Vdia,i "b v poi,i H Jj (4) 

Ti 1 dni^ 1 d — V</> x n 

]e\Bn' l ~dT ev ~ BT^di ^' + B 
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The equation for the velocity of the electrons is 

V e = Vj_,e +V|| >e 



V± e = V E + V dla . e (5) 

- V(A x n 11 



B n rripfl 



(n x Vp) 



-V<ft x n | T e 1 dn ^ 

B —\e\B n dr &y W 

We assume neutrality n e = rij = n and introduce the expressions of the 
velocities in the continuity equations for ions and for electrons. 
The ion continuity equation is 

— + V-(nv ±>< ) = 
The electron continuity equation is 

^ + V_L- (nV_L, e ) + V|| (nU|| e ) = 

The equations are substracted , to obtain 

-V|| (nu|| e ) 

n—'V^cf) 



BQ., \ dt 
+V ± • (nv di0)i ) - V ± • (nv diQie ) 




Ti ( T e 1 dn \ dn 



From the last term in the left we get: 
■^ej • Vdia,i + ( V±n) • v dw . 



T e \ |e| £? n rfx / 9y 
Including the similar term for the electrons, we obtain 

-V|| (n«|| e ) 

1 „ ( d ^ 
-V_l • n—V_ 



Bfli V dt 
Tj ( T e 1 (inn \ <9n / T„ 1 rfrin \ 9n 



T e \|e|_Bno rfa; / dy \\e\Bno dx J dy 

= 
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or: 



- l + =f 



T,; \ f T e 1 dn \ dn 1 



T e J \\e\ B no dx J dy Bfli 



From the continuity equation for electrons 
d 



we obtain 



— + V± • (nvi) + V|| • (nv||) = 



dn 1 90 dn „ <9n 

"^7 + ~5 a 3 h +Vdia.e ' V _L?l + Vq — 

at B ay dx ay 



(7) 
(8) 
(9) 



-n (V|| • v||) + «||V 







where the seed poloidal velocity is Vo (x) — — d< ^^ e y . The parallel moment- 
tum balance gives the parallel electron velocity 



v He = -^ V ll (|e|^ + T e ln(n/n )) 



(10) 



In the absence of friction (a — > oo) and of particle drifts the electron response 
is adiabatic 



n 
no 



and the potential is determined from Eq.(^). To develop separately the ion- 
polarization drift term, we introduce the notation: 



W 



1 



Bn 
i 

bU 



(mi - ( • v/, • V ) V 



■V 

-V_l- [(n + n)I 



where we make explicit the electric potential cj> associated to the initial plasma 
poloidal rotation, Voe y . 



d_ 

dt 



(V e y + v) • V H V_l (<£„ + 0) 



We have 



dt 



'at 



G 



or 



with the relations 



— (e x BV a ) + 



i M 



+ 



dy 



(v ± 4> + v±0) 



at ay 

+ ^(v^)+(v-V ± )W e x + 

+%|(v^) + (^v 1 )(v 1 ^) 



V_l0 = B(vxn) 

= + (-Bv x )e y 



~ ld< t> A~ 
v * = Bte aDdv * 

After very simple calculations we obtain: 

1 



1 90 
5^ 



. ~ 9Vq 9Vq 
at ax ay 



dx 



+ v. 



OVy 

V dy 



and 



It will be useful to calculate the derivaties of these quantities 



}_dh = d dV Q | dv v dV | _ d 2 V | dVp 8Vq | T , d 2 V 
B dx dx dt dx dy v dydx dx dy dxdy 

dv v dV „ d 2 V d dv v 

_| » _|_ y _ _| Z. _|_ 

9i 9a; 21 fe 2 dt dx 

dV OVy d 2 Vy 

Vvi — 



dx dy dxdy 
dv x dv x „ d 2 v y 

H ^ ^ ~t~ T> o~ ~t~ 

os dy dx z 
dv v dv v „ 9 2 w„ 



and 



9a; <9y <9a;<9y 



1 dl y _ dV dv x 
B dy dy dy 



d dv x d 2 v % 



dt dy dy 2 

dv x dv x „ d 2 i 

- v. 



dy dx x dxdy 
dv y dv x „ d 2 -^ 



9y cfy w 9y 2 
The quantity denoted by W takes the form 

w = mvte) Ix + m no {^) + 

1 /0n\ 7 . 1 



(14) 



BSli \dx J Bili \ dx 
i /av\ J_ f&n\ 
+ Bn t n ° \ dy ) + Bn t \dy)* v ' " V 0» 

2.1 The mode evolution in a fixed plasma rotation profile 



We will assume that the mode evolves initially without perturbing the equilib- 
rium profiles, in particular the seed poloidal rotation. This allows us to simplify 
the expressions above, taking: 




= 
= 



8 



Then the first lines in the formulas Eqs.Q), @, Q are zero. Let us consider 
in the expression of W the part Wq which does not contain the fluctuating 
density n. Writting 



we have 
W 











i 









W — Wq + W 

dv„ _ 9tL 



9V , T/ 
aa; ay 



(9? 



9v„ 



and 



W = 



1 



9n\ 



BQi \dx J 



dx 



Ba, 



dx 



On 
dy 



dy dt 



BQi 



dly 

dy 



Replacing the perturbed velocity by the perturbed potential, writting all terms 
and summing, we get: 



Wo 



d_ 

ot 

l 



V 



■no 



0_ 

dy J B dx 



1 

B 



-Vj_0 x n 
B 



1 90 \ dV 
B dy I dx 



-V^</> x n 
B 



90 
dx 



1 d(j> d 2 V 
B dy dx 2 

Collecting all what we have at this moment the ion continuity equation (0) 
becomes: 



V,| (nv\\ e ) + 



T e 1 dno \ dn 



1 



5ft 



■n 



e\B no dx ) dy 
d 



Vo 



dy 



Y72 1 \r" d ^ 



-Vj_0 x n 



1 



da; 
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dt 



dy 



9v„ 



v: 

dv v 



dV 



tt: + Vfo — v v + v x — h u,,— h w^— — 



9a; 



9y a 9a; 



+ {terms of the first line in the expressions of I x and derivatives} 



= 



In the above equation (which is exact) we shall make the following approxi- 
mations: 



9 



• neglect the term containing the parallel electron velocity, assuming infinite 
electric conductivity; 

• neglect the term which contains dno/dx since it is in the ratio k : 1/L n 
with the other terms, and we consider 

kL n < 1 

• neglect W ; these terms are in the ratio n/no with the terms which are 
retained; 

• neglect of the terms in the first lines of the expressions for I x and dl x /dx, 
dly/dy. (These are the terms in the curly brakets, the last line). As 
explained above, we assume that the mode evolves in a background of 
fixed rotation profile, Vq (x). 

The resulting equation is 

i} + T^) (|e|_Bn dx ) dy (n ) 

+ U °d~y) ±<P ~ °d7j \ B—' v± ) v ^ 

= 

and, replace the adiabatic form of the density perturbation 

n \e\ </> 
n T e 

We define 

/ TA ( T e 1 dn \ lei / E\ 1 

and obtain 

80 (d , v d\ 2 ~ „dt / -Vi^xn \ 2 ~ 
which is the barotropic equation. 

2.2 Nondimensional form of the equation 

We consider that the ion mode extends in the spatial (x) direction over a length 
L. A typical value for the sheared poloidal rotation is noted Uq- We make the 
replacements 

V -» Ly 
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V - UU 

such that from now on y, t, <f) and U arc nondimensional quantities. We also 
change the radial coordinate into a dimensionless variable 

x — ► La; 

and rewrite the equation 

TA 1 i 2 



dx 2 dy 





The coefficients are 

/ T\ 1 L 2 
= Oj I 1 + — % 7 I — — — non-dimensional 

\ TeJLn Uo 

e = — ?4 ^" non-dimensional 
f7 |e| B L 

For an order of magnitude, e is the ratio of the diamagnetic electron velocity 
to the rotation velocity Uo multiplied by the ratio of the density gradient length 
to the length of the spatial domain. This quantity, e is in general smaller than 
unity. 

The quantity is the ratio of the ion cyclotron frequency to the inverse 
of the time required to cross the spatial domain with the typical flow velocity. 
Since the later (Uq/L) involves macroscopic quantities this ratio can be large. 
It is multiplied by the ratio of the spatial length to the density gradient length 
(these quantities can be comparable and the ratio not too different of unity) . 

We change the notations eliminating the primes. The equation becomes 

^^(l^|)^-Sf^ + " ( -^ xa, - VJV1 *- 

This is the barotropic equation, known in the physics of the atmosphere. 
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3 Multiple time and space scale analysis of the 
ion mode equation reduced to the barotropic 
form 

The quasigeostrophic barotropic atmospheric model leads to the barotropic 
equation (see Horton M where the quasigeostrophic approximation is explained 
in relation with the Ertel's theorem). This equation has been studied in the 
atmosphere science by means of the multiple space and time scales method. Our 
aim is similar, i.e. to obtain an equation for the envelope of the fluctuating held 
4> in order to study the possible poloidally nonuniform profiles of the turbulence 
averaged amplitude. The connection between the original equation and the final 
slow-time and large-spatial scales equations is not simple : much numerical work 
is required in order to connect the input physical parameters of the plasma 
with the formal coefficients in the final equations. The detailed analytical work 
is presented in Ref(|^]) and the necessary formulas are reproduced here for 
convenience. To simplify the expressions the following notation is introduced 
[4>, V 2 </>] = [(— x n) • Vj_] V^_</> and the equation can be written 

The analysis on multiple scales starts by introducing the variables on scales 
separated by the parameter e: 

Ti = et, T 2 = e 2 t 
Yi = ey, Y 2 = e 2 y 

This gives for the derivatives 

d d d 2 d 
di ~* dt +£ Wi +£ df 2 

d d d 2 d 
dy dy dYi dY 2 

The solution is adopted to be of the form 

We denote the linear part of the operator in the equation by 
( d ,J\ 2 / d 2 U(x)\ d 
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(15) 



(16) 
(17) 



Substituting the Eq.(|T7|) and the forms of the derivation operators Eqs.(|l5|, |l6| ) 
in the equation of motion we obtain from the equality between the coefficients 
of the powers of the variable e: 



L(j> 



(i) 







(18) 



(2) 



_0 



(i) 



d 2 ^ 



d 



d 



dt dy J dydYi 



(19) 



L(f> 



(3) 



d_ 

dt 



dY 2/ 

d \ d 2 (j> {1) 



U 



dy J dydY 2 
9V« 



d_ d 

dt + dy) dY 2 



d 
dTi 
d 



U- 



d \ d 2 <t> {1) 



dYi) dydYi 
dTi dYj v 



d 



dY 2 



d \ d 2 (h {2) 



dt dy J dydY] 



dYi 



dYi dx dx 8Y\ 

2 ( d<t> {1) d^ {1) _ d^ a 3 (1) \ 

I dy dydxdYi dx dy 2 dY\ J 



-,(2) V7 2 ^(!) 



(20) 



The solution is represented on the space and time scales. 
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3.1 The enveloppe equation 

The solution of the linear part of the operator, Eg. is adopted as the super- 
position of two propagating waves 

4> {1) = Ax(T u T 2) Y 1 ,Y 2 )(p 1 (x)exp(ikiy-iuj 1 t) (21) 
+A 2 (Ti,T 2 ,Yi,Y 2 )ip 2 (x) exp (ik 2 y - iu 2 t) 
+c.c. 

The two functions </?„, n = 1,2, are solutions of the equation 

-d^-\ K ~ u^r)^^ (22) 

with 

Vn (0) = ip n (L) = 

In this formula, (0, L) are the limits of the spatial domain in the minor radius 
direction where the flow U (x) ^ ; c„ = u> n /k n are the phase velocities. To 
find the solutions of this equation we have to use numerical methods. This is 
a Schrodinger-type equation where the potential depends on the energy, (since 
the phase velocity c n depends on the wavenumber k n ). We must use numerical 
methods suitable for Sturm-Liouville equations and in general the solution exists 
for only particular values of the parameters. We choose to fix the frequencies u> n 
and consider k n as the parameter to be determined. The problem is complicated 
by the possible occurence of singularities in the potential due to resonances 
between the phase velocity and the fixed zonal flow. We will avoid any resonance 
since the physical content of the processes consisting of the direct energy transfer 
between the fluctuations and the flow represents a different physical problem 
and requires separate consideration. In the present work we examine the effect of 
the turbulence on the poloidal rotation as being due to the nonuniform diffusion 
rates and Stringer mechanism. 

The parameters in the Eqs.(|22]) are f3, e and the sheared flow velocity U(x). 
We take for U (x) a symmetric Gaussian-like profile shifted in amplitude and 
retain its maximum, half- width and asymptotic value (i.e. at x — or L) as 
parameters. The two frequencies have fixed values during the eigenvalue search 
but they must be considered free parameters as well. This builds up a large 
space of parameters which should be sampled. For some point in this space the 
eigenvalue problems Eqs.(^||) are solved and the phase velocities compared with 
the profile of U (x) . The presence of resonances renders the set of parameters 
useless. Since (as will result from the formuals below) there are also other 
possible resonances, the determination of a useful solution (kx,k 2 ) is difficult 
and requires many trials. 

The numerical methods that have been used in solving Eqs.(j2^) belong to two 
classes: boundary value integration methods and respectively shooting meth- 
ods. When the sign of the "potential" in the equation is negative everywhere on 



14 



(0, L) , which occurs for short wavelengths and smooth fluid velocity profile the 
boundary conditions cannot be fulfilled except for asymtotically, i.e. assuming 
arbitrary small nonzero values at these limits, since the solution decays expo- 
nentially. However most of the boundary value integrators we have tried gave 
amplitudes of magnitudes comparable to the boundary value, which essentially 
are homotopic to the trivial vanishing solution (the equation is homogeneous) 
and render useless this approach. We have to use a shooting method accepting 
the difficulty that the initial value of the parameter to be determined (eigen- 
value) must be placed not far from the final value. We have used the specialized 
package SLEIGN2 interactively. For the investigation of ranges of parameters 
most of the calculations have been performed using the routine D02KEF of 
NAG. For identical conditions the two numerical codes gave the same results 
within the accepted tolerance. 

Using the routine D02KEF we divide the integration range (0, L) into four 
intervals corresponding roughly to a different behaviour of the potential. This 
should simplify the task of integrator. A single matching point is assumed, 
usually at the centre. The tolerence is 10~ 8 which however does not exclude 
fake convergence in case of an initialisation of the eigenvalue very far from the 
true solution.. 

We notice from Eqs.(p2[) that the potential can easily be dominated by high 
shear U" especially for smaller wavenumbers k n . Actually we are more interested 
to investigate regimes where only mild effects of the seed sheared rotation arise 
during the process of self-modulation since, as explained before, the possible 
evolution to the nonuniformity in the poloidal direction will be the source of the 
torque on the plasma. Avoiding regimes of [/"-dominated potentials we assume 
low shear and take the Gaussian profile with U max — E/min = 2 x 10 4 (m/s) 
represented in Figure (|l|). The following set of physical parameters has been 
used: minor radius a = 1 m, magnetic field Bt = 3.5 T, electron temperature 
T e = 1 KeV, ion temperature T { = 1 KeV, density n = 10 20 m~ 3 , L n = 10 m. 
The parameters in the barotropic equation are then: (3 = 33.526 and e — 
0.28571. The two frequencies are u x = 85 x 10 2 s" 1 and lo 2 = 85 x 10 4 s -1 . It 
is assumed that the space region of significant magnitude of U (x) is L = 0.1 m. 
The following two eigenvalues are obtained, normalised to L _1 : fci = 0.5525 
and &2 = 3.2029. The potentials in the Schrodinger-like equations and the 
eigenfunctions are shown in Figure (^|). 

The two eigenfunctions (fi 2 have an important role in the following steps 
of the calculation. Substituting the Eq.(21) in the barotropic equation and 
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expanding the operators, we obtain: 

2 



L0 (2) = - Gi n exp (ik n y - iu n t) 



where 



We have fci 



n=l 
2 



ig2nA 2 n exp [£2 (k n y - u) n t)} 



n=l 



-g 3 AxA 2 exp [i (k 12 y - Wi 2 *)] 
-ig±A Y A* 2 exp [i (a 12 y - o 12 i)\ 
+c.c. 



ku = ki + k 2 

£J 12 = Wi + U! 2 



3.755, w = 8.585. 



0112 
C12 



fci - fc 2 



or: ai2 

Gin — 



2.650, (Ti 2 

dx 2 
U 



.415. 



dAn 



dx 2 



(23) 



-k 2 n <p n )+(f3-U")<p n -2k n (U-c n ) 



dAn 

m 



g 2 =k (ip — - —\ d2lfn 
n \ n dx dx J dx 2 

These functions can be obtained after calculating the eigenfunctions tp 1 2 and 
are represented in Figure (|^). 
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The function g 2] 




0.2 0.4 0.6 0.8 1 
The function g 





0.2 0.4 0.6 0.E 
The function g 




0.2 0.4 0.6 



Figure 3: Graphs of the functions <?'s. 



Various solutions to the equation ( p3| ) are possible, taking the space-time 
dependence similar to one of the terms which compose the right hand side. 
The second term suggests solutions of the type 

02n = W 2n Al exp [i2 (k n y - oj n t)] + c.c. 
The third term suggets solutions of the type 

(2) 

03 = W 3 A 1 A 2 exp [i (k 12 y - iv 12 t)] + c.c. 
The fourth term gives 

(2) 

4>\ = W±AxA* 2 exp [i (ai 2 y - <Ti 2 t)} + c.c. 
In the formulas above the functions W 2n (x) are solutions of the equations: 

92n 



d 2 W 2n f Al2 0-U" 



dx 2 



d 2 W 3 



dx 2 



d 2 Wj 
dx 2 



k 2 

ft 12 



W' 



2/i 



U - C r , 

with W 2n (0) 



W 3 = 



(£ - U") 
Uk\ 2 — ^12 

with W 3 {0) 



an (J3-U") 



l \2 



Uot\ 2 — <7i2 

with W 4 (0) 



2 (Uk n - UJ n ) 

W 2n (L) = 

93 

Uki2 - UJ\2 

W 3 (L) = 

U0L\ 2 — 0\ 2 

Wi (L) = 



18 




Figure 4: Graphs of the potentials occuring in the equations for the functions 
w's. 



These equations are integrated numerically with a boundary value integrator, a 
finite difference method. The results are shown in Figures (jij), (^) and (|^). In the 
following figures are explained also the functions W\ n which will be introduced 
in Eq.@. 

The first term gives a solution of the following form 



A 2 ) ( 2 ) i-i ■ +\ , 

4>\n = 'Pin ex P \ lk nV ~ luJ nt) + C.C. 



dx 2 



S 2 ) 
'in 

p-U" 



where the function </?^ n satisfies the equation 



u ■ 



(2) 
Pin 



U 



dx 2 



-2k 2 n (U-c n )] 



dx 2 




\ dA n 
1 8T X 


k l<Pn^ 


+ (13- 


U»)<p 










\ 





This equation serves to generate a condition on the two amplitudes A\_2- 
When the left hand side of the equation is multiplied by the function ip n 
and integrated between the limits on the x domain, it gives zero, due to the 
boundary conditions. The same must then be true for the right hand side. This 
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Figure 6: Graphs of the functions w's. 
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gives a "solubility condition", on the space time scales slower at order 1: 



dY 1 



where 



yn 



knVndx/ 







13 -U" 



o (u- c n y 



-ip 2 n dx 



Expressed in units of length L and units of time to = 10 4 , the group velocities 
are : c gl = 0.30328 and c g2 = 3.9908. 

Then it is clear that the amplitude A n propagates at the speed c gn : 

A n = A n {Y x -c gn T x ,Y 2 ,T 2 ) 

Now we introduce a new function 



(2) 
fin 



The equation satisfied by W\ n is 

d 2 W ln f l2 P-U" 



dx 2 



i 



(c„ 



-gn ) 



iW ln 

Wi n 
dx 2 



dA n 
dY x 



(24) 



with the bounday conditions 

W ln (0) = Wm (L) = 
At this order (2) there are the following solutions 

^(2) _ 



Kfn ~ 2k n (U ~ Cn) ip n 



4? 



,(2) 
'22 



(25) 



0f 



+*(ar,Ti,yi) 



The last function will be determined by the condition of solubility at the next 
order on the space-time scales (O (e 2 ))- To obtain the evolution equations for 
$ and A n , we use the expressions for Eq.(pl|) and 4>^ 2 \ Eq.(p5|) in the right 
hand side of the equation (20), and all inhomogeneous terms are obtained. In 
additions there are terms which are not dependent on y and t, and these terms 
should be zero. This leads to the solubility condition relating the function $ 
(correction to the mean flow due to the presence of the wave packets) to the 
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amplitude A n . 



m +u w 1 )^ + ^- u) w 1 



(26) 



2 



J 1 da; 2 

d 



d<y5„\ / dV 



dx dx 



dx 2 



n ^ n dx 



This relation will be used later. 

Returning to the equation (^0|), we examine the content of the inhomoge- 
neous part (the right hand side). There are terms which have a space-time 
dependence 



exp [ik n y — iu n t] 



(27) 



and in order to have solutions of this type, we need to impose solubility condi- 
tions, otherwise there will arise secular terms at the order To find these 
restrictions (solubility conditions) we multiply the equation with the function 
conjugated to ( |27| ) 

exp [-ik n y + iuj n t] 

and integrate time on the interval (0,2%/u) n ) and space: y on the interval 
(0,27r/fc n ) and x on the interval (0, L). On the left hand side we obtain zero, 
so this is what we must obtain also from the right hand side. The conditions 
which results are: 



d 







+ c„i—— A\ - iax 



dY 2 



i(p 1 |A 1 | 2 + 7l2 |^ 2 | 2 + A 1 )A 1 = 



-^ + c g 2^p)A 2 -i a2 ^0-i(p 2 \A 2 \ 2 + l21 \A 1 \ 2 + \ 2 ) A, =n 



dY 2 



The notations which have been introduced are: a n = If 2 -, p n — -jfS 7i 2 = nf j 
721 = TT7' A ™ = iff and 



n„ = - 



U -c, 



-fndx , I n = 



U -c, 



-fnodx 



U-c, 



<Pi 



o U -a 



f 12 dx 



hx = 



<P 2 

o U - c 2 



f 21 dx , I n3 = 



U-C, 



-fnzdx 
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The equations satisfied by all these functions are given given in Ref.([H). 



fn 



dx 2 



-U 



d 2 fn 

dx 2 



-k 2 n <p n )-(f3-U")<p n + 2k 2 (U-Cn) 



fnO = -2k n (U - C gn )(f n - k n (U - C n ) 



(U - c gn ) 



yn) 

, d 2 W 1 , 
' ^ dy 2 



-klWv 



+2k 2 n (U - c n ) W ln - (/3 - U") W ln 



, w d dW 2n 
2k n W 2n — + k n—, — 
dx dx 



dx 2 



k l<Pn 



•2n 



dx 2 



2n 
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, u . d dW 3 
ki 2 W 3 — + k 2 —— 
dx dx 

dx dx 
d dif 2 

k 2H>2-, "12-^— 

dx dx 

w d u dWi 
ai 2 Wi- k 2 —— 

dx dx 



k 12 W 3 — + k 1 
dx 



dW 3 
dx 



d dip-, 

kiVx— + a i2-^— 
dx dx 

t d . , dip 2 
kiip 1 — + ki 2 —— 
dx dx 

„, d , dW 4 
dx dx 



d 2 ip 2 
dx 2 



dx 2 
d 2 W 4 
dx 2 

d 2 y> 2 
dx 2 

d 2 <Pi 
dx 2 
?W< 
dx 2 

fw 3 

dx 2 

d^fi 
dx 2 



k\ip 2 



k 2 W 3 



a{ 2 W A 

- k 2 tp 2 

- k 2 (f! 
- a\ 2 Wi 



k\ 2 W 3 



fn3 



kn 1 



d 2 f n 
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Figure 7: Graphs of the first group of functions /. 




24 



We use a boundary value integrator for the inhomogeneous equations and (the 
NAG routine D02JAF). The results are given in the Figures (0) and 

After numerical integrations the follwing values are obtained for the con- 
stants: a-i = 0.39976, a 2 = 0.45246, p l = -0.10182, p 2 = 1.8115, 7 12 = 
0.091474, 721 = -5.3768. 

The following transformation (Jeffrey and Kawahara) is performed in Ref.( 



T = T 2 

X = - £ {Y 2 -c gl T 2 ) = Y 1 -c gl T 1 



Then the equation of connection between $ and A n becomes 

02. 

dx 

dx 



_d_ 

dY 

2 

E 



n=l 



/ d2 (vsr d( ?n dw m \ A , d<p n 



- V 



dl Pn\ ( d2 Vn 



1 dx dx I V dx 2 
The solution of the equation ( pq ) is considered to be of the form 



(28) 



$ = H n (*) \A n \ 



n=l 



and the equation for H n (x) results: 



(U-c gl )^ + (f3-U")H n 



, d 2 ( dip,, dW lr 



(fa: 



e? d</?„ \ f d 2 



t ^ n dx dx 



dx 2 



(29) 



The change of variables is performed in the equations for the amplitudes A r 
and, again, the Eq.(|29|) is used: 



dA 1 d 2 A l 



dT 



dY 2 



U \A 



U +isi2\A 2 \ 2 ) Ai =0 



, d . d \ . d 2 A 2 
l {df- 5 dY) A2 + a2 ^ 



(T2\A 2 \ 2 + V 2 l\A 1 \ 2 ) A 2 =0 



where 



5=~ (Cgl ~ C g2 ) 
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Potential for H. 



The function H 




Figure 9: Graphs of the functions H. 



n n Jo u — c n 



d 3 H n 
" dx 3 



dx 2 



dHn 

dx 



dx 



V12 = 7l2 



1 f L fciy 

niio u- 



ci 



<Pi- 



d 3 H 2 
dx 3 



dx 2 



dx 



dx 



V<1\ 



721 



k 2 ^ 2 

n U 



C2 



<^2' 



cix 3 



dx 2 



The following transformation is convenient: 



A 2 = B cxp 



2a 2 



y 



This leads to the following equations: 



.8 At 



d 2 A l 
dY 2 



dx 



dx 



(aMA 2 + v 12 \B\ 2 ^ A 1= 



(30) 



§ + a 2 ^ + (a 2 \Bf + , 21 \A 1 f)B 







(31) 



The constants are: a\ = —1.6324, a 2 — 1.9331, v\ 2 = —0.53387, v 2 \ = 
—2.9014. The constants ai j2 are called the "dispersion", 171,2 are called the 
nonlinearity ("Landau") constants ; 1/12,21 represent the coupling of the two 
amplitudes. 
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Numerical experience shows that the dispersion constants are sensitive to 
the precision of the integration scheme. We use for D02KEF the value for 
tolerence 10 -8 . 

The system of coupled Cubic Nonlinear Schrodinger Equations has been inte- 
grated numerically in relation with many applications p0| , p"l| ] . One important 
conclusion is that, in the case of a single equation, the sign of the product of 
the dispersion and of the nonlinearity parameters establishes the possibility of 
a soliton formation. 

3.2 Numerical search for classes of admissible envelope 
solutions 

In principle an interactive search for eigenvalues (either with SLEIGN or D02KEF) 
is able to obtain the eigenvalues k\ , eventually after series of less-inspired 
initializations. However in many situations resonances occur and the sources 
and/or the potentials in the differential equations for the functions W n , f n be- 
come discontinuous. Since these eigenvalues cannot be used we have to devise 
an automatic search over the space of parameters with checks of continuity of 
the functions and automatic discarding of wrong results. This is accomplished 
by the routine P01ABF of NAG. 

The search reveals that the short and respectively long radial wavelength 
part of the ion wave spectrum have different behavior under the modulation 
instability. The most affected is the long wavelength part (small k x ) with u> 
of the order of the ion-diamagnetic frequency. This corresponds to the domi- 
nant part of the ion spectrum, because, as can be seen in numerous numerical 
simulations, the dominant structures are elongated in the radial direction, with 
various degrees of tilting compared to the radius. It is also important to note 
that finding the eigenvalues in the scan of the parameter space weas possible in 
the condition that the eigenfunctions ip 1 2 (x) have low number of oscillations in 
the radial direction (one or two) . Again this corresponds to the important part 
of the ITG turbulence spectrum. 

The variation of the parameters of the system of NSE with the physical 
parameters of the problem can be studied by the automatic search of eigenvalues, 
as described before. Below are plotted four graphs showing this dependence. 

4 Exact solutions of the Nonlinear Schrodinger 
Equation and the nonlinear stability problem 

4.1 Introduction 

The Nonlinear Schrodinger Equation can be solved exactly on an infinite spatial 
domain using the Inverse Scattering Transform. The first step is to represent the 
nonlinear equation as a condition of compatibility of two linear equations. This 
is achieved by finding a pair of linear operators (Lax operators) and noticing 
that for the first of them the eigenfunction problem is a Schrodinger equation 
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Figure 10: Parameter dependence of the coefficient a\. 





Figure 12: Parameter dependence of the coefficient <j\. 
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leading to a quantum scattering problem. The unknown function of the NSE is 
the potential of the Schrodinger equation. Using the initial condition of the NSE 
the scattering data are determined at t = 0. The scattering data have simple 
time dependence and so they can be evloved in time to the desired moment. Re- 
turning from this new set of scattering data (to the potential that has produced 
it) is achieved by solving an integral equation (Gelfand-Levitan-Marchenko) and 
provides the solution of the NSE. In general this solution consists of solitons and 
"radiation" . 

The Inverse Scattering Transform on periodic domain [|l4| , 1 13 , Jll| is a more 



powerful procedure since it reveals and takes advantage of the deep geometric 
and topological nature of the problem. The admissible solutions of the Lax 
eigenvalue problem can now only be periodic functions. Since after one period 
the change of any solution can only be linear (via the monodromy matrix) 
the periodic or antiperiodic functions must be found as eigenfunctions of this 
monodromy operator (matrix), under the constrain that the eigenvalues are 
complex of unit absolute value. This singles out a set of complex values of 
the spectral parameter A (the formal eigenvalue in the Lax equation), the main 
spectrum. From this set, the particular values that make the two eigenfunctions 
to coincide are called non-degenerate. Here the squared Wronskian (which is 
a space and time invariant) has zeros on the complex A plane. Then the non- 
degenerate A's are singular points for the Wronskian. Removing the square-root 
indeterminancy one has to connect pairs of non-degenerate eigenvalues by cuts 
in the A plane. The Wronskian becomes a hyperelliptic Riemann surface. The 
evolution of the unknown NSE solutions (called hyperelliptic functions /x - (x,t)) 
on this surface is as complicated as the original NSE equation. However, via 
the Abel map the Riemann surface is mapped onto a torus and on this torus the 
motion magically becomes linear. After obtaining the space-time dependence 
on the torus we need to come back to the original framework. This is called the 
Jacobi inversion and can be done exactly in terms of Riemann theta functions, 
giving at the end the exact solution of the NSE. 

Not only the geometric approach is more clear but it also allows a treatment 
of the stability problem for the solutions of the NSE since it allows to trace the 
changes of the main spectrum after a perturbation of the initial condition. The 
following subsections provide a more detailed discussion of the 1ST of a periodic 
domain. The next section will discuss the stability. This information is available 
from the abundant literature on the 1ST and is only mentioned here in order 
to understand the mechanism governing the stability of the envelope of the ion 
wave turbulence. For this reason we will focus on the determination of the main 
spectrum and its role in the construction of the solution. The effective steps to 
be undertaken to obtain the solution will only be briefly described. We strongly 
recommend the lecture of Ref.(|l^]) on which this presentation is based. 

4.2 The Lax operators and the main spectrum 

The 1ST method starts by introducing a pair of linear operators, called the Lax 
pair (see [O) which allow to express the nonlinear equation as a compatibil- 



30 



ity condition for a system of two linear equations. For the cubic Schrodinger 
equation 



.du d 2 



i- 



U , nl |2 



dt dx 2 
the Lax pair of operators is 



2\u\u = Q (32) 



is( !t rt ) (33) 

1 — u [x, t) —to x 1 



_ , i\u\ 2 -2i\ 2 -u x + 2i\u , ,. . 

1 i£ + 2iAu* -*|t*| a + 2*A a 1 



The action of the operators on a two-component (column) wave function 
is given by the equations 



>2 



L<p = Xcj) (35) 

and the condition of compatibility of these two equations <f> xt — <f> tx is precisely 
the cubic Nonlinear Schrodinger Equation. 

4.3 Solving the eigenvalue equation for the operator L 

As in the "infinite-domain" 1ST, we have to solve the eigenvalue equation 
Eq.(p5|) using the initial condition for the function u(x,t), u(x,0). How- 
ever, for the particular case of periodic solutions of the NSE, we must have 
u {x + d, 0) = U (x, 0) where d is the length of the spatial period (i.e. actually L 
in the previous notation; however we use d in this and the next section, keeping 
L for the Lax operator). Fixing an arbitrary base point x — xo one considers 



the two independent solutions of the equation (35) which take the following 
"initial" values at x — xo: 

= ( J ) and ^(*o) = ( ? ) (37) 
The matrix $ (a;, xq; X) of solutions is given by 

$ (x,x ; X) = ( ^'^1 ) ( 38 ) 

\ 02 [%,%o;X) 2 (x,x ;X) J 

i.e. this solution satisfies the equation 

L$ = A$ and <$> (x , x ; X) = ( q i ) ( 39 ) 
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It can be seen that the Wronskian of f 0, 0J is the determinant of the matrix 
of fundamental solutions 



W 



(0,0) =det($) 



It can be shown that = i.e. the Wronskian is constant, which gives 

det($(x)) = det($(a;o)) = 1 

Using the initial condition for <E> one can solve the equation (^) and obtain 
$ (x, xq; A); in particular we can obtain it at xq + d, i.e. after a period : 

$ (xq + djX0 . x)= ( <h to + d,xa; A) £ (xo + j^o;A) \ 
\ 02 to + ">, x ; A) 02 to + a, x ; A) J 

This is called the monodromy matrix and is noted 

A/(x ;A) = <5>(x + d,x ;\) (41) 

The monodromy matrix is the matrix of the fundamental solutions calculated 
for one spatial period. In general the monodromy matrix is the element of the 
monodromy group associated to a loop and to its homotopy class, for a marked 
point on a manifold (here the circle S 1 = [0, d). The linear monodromy group 
acts by replacing the elements of the vector column by linear combinations 
of the initial ones. For changes of the marked point xq the matrix has two 
invariants : the trace and the determinant: 

[TrM] {x ;X) = [TrM] (A) = A (A) 



det M = 1 

The function A (A) is called the discriminant and is independent of Xq- 

Looking for Bloch (Floquet) solutions The values the discriminant A (A) 
takes on the complex A-plane control the monodromy and by consequence select 
the values of A for which admissible eigenfunction of the Lax problem exist. 
In other words A governs the spectral properties of the operator L. To find 
under what conditions periodic solutions are possible, we construct the Bloch 
(or Floquet) solutions of the equation L<j) = A0. The Bloch function has the 
property 

i(> (x + d; A) = e l ^ A V (x; A) (42) 

where p (A) is the Floquet exponent. Like any other solution of the Lax 
eigenproblem tp can be expressed as a linear combination of the two fundamental 
solutions <j) and 

■0 (x; A) = A<p (x; A) + 50 (x; A) 
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For the particular point xq we have, taking into account the "initial" conditions 

ip (x ; A) = Acp (x ; A) + Bcp (x ; A) = 

After one period d the function ip is linearly modified by the monodromy matrix 
and, according to Eq.(f42]) is multiplied by a complex number of unit absolute 
value. We write m (A) for this number and look for those A's where this is 
complex of modulus one. The regions of the A plane where to (A) has modulus 
different of one correspond to unstable functions ip. We write 

ip(x +d; A) = to (A) ip(x ; A) (43) 




M U) =m(A) U) (44) 

The Bloch functions are invariant directions for the monodromy operator M 
acting in the base formed by the fundamental solutions of the Lax eigenproblcm. 
The monodromy operator preserves these directions and multiplies the Bloch 
vector with m (A) = exp (ip (A)). 

det (M -m)=m 2 - (Tr M)m + dct M = m 2 - A (A) to + 1 = 

This gives 

A(A)±(A 2 (A) -4) 1/2 
m (A) = y — '- — 

In general A (A) is an analytic function on the plane of the complex vari- 
able A. The equation Im [A (A)] = is a single relation with two unknowns, 
(Re A, Im A) and calculating one of them leaves the dependence on the other as 
a free parameter. This gives a curve in the plane and the real-A axis is such a 
curve: for Im A = we have Im A (A) = 0. 

We look at the variation of A (A) on the complex plane to find the effect on 
m ± (A). 

• Where the Floquet multiplier is a complex number of magnitude unity 
the functions are only effected by a phase factor after one period. The 
Bloch functions (eigenfunctions of L) are stable under translations with 
d. The real A axis is a region of stability. The region where A is real and 

A (A) < 4 

is the band of stability, the modulus of to being 1. 

• When A is such that the discriminant (A (A) = Tr (M) ) is equal to 4, 
the eigenvalues of the monodromy matrix, to (A), are ±1, and the Bloch 
functions are periodic or antiperiodic. This is the Main Spectrum in A 
and is noted {A.,}. The main spectrum consists of a set of discrete values 
A,-. 
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• for those A which gives \m (A)| ^ 1 the Bloch eigenfunctions are unstable. 

Consider the case where A (A) is complex. For all values of the parameter A 
on the complex plane, [with the exception of the Main Spectrum where A (A) = 
±2 ], the eigenvalues of the monodromy matrix are distinct which means that 
the eigenfunctions Bloch arc distinct and Independent (The Wronskian is 
different of zero) . 



However we can only be intersted in points of the main spectrum, as they 
can provide admissible (periodic) eigenfunctions of the Lax operator. (1) There 
are points Xj in the main spectrum where the two eigenfunctions Bloch are 
distinct and independent: these points are called degenerate (with the meaning 
that for two independent eigenfunctions there is only one eigenvalue). (2) There 
are points Xj where the two eigenfunctions Bloch are NOT independent: these 
points are called nondegenerate. For such values of A the Wronskian is zero. 

The important class of initial conditions u(x,t — 0) with the property that 
there is only a finite number of non-degenerate eigenvalues Xj is called finite- 
band potential. 

4.4 The "squared" eigenfunctions 

Starting from the initial condition u (x, 0) and leting it evolve in time according 

to the NSE, u(x,0) N -^ E u(x,t) 1 the discriminant remains invariant. All the 
spectral structure obtained from the discriminant A (A), i.e. the main spectrum, 
the stability bands, are invariant. The eigenvalues of the monodromy matrix Xj 
are invariant. 

4.4.1 The two Bloch functions for the operator L 

Consider the two Bloch eigenfunctions of the operator L: 



(x + d) = m + ip + (x) 
(x + d) = m~ip~ (x) 




and define the following functions 



f(x,t; A) = 



+ 



(45) 



g{x,t;X) = ViVi 



(46) 



h {x, t; A) 



(47) 
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These functions have the property that they are periodic, / (x + d, t; A) = 
f (x,t;\)foT all values of A. Since / is composed of ip + and ip~ and these 
are Bloch solutions for the eigenvalue equation Lip = Xip we have 

^L=u* g -uh (48) 
= -2i\g - 2uf (49) 

dh 

— = 2iXh + 2u*f (50) 
ox 

and similarly we write the time derivatives taking account of the compatibility 
condition which involves the operator A. 

df 

J - - {u* x + 2iXu*) g + i (-u x + 2iXu) h (51) 



dt 



dt 



2i(\u\ 2 g + 2i(-u x + 2i\u) f (52) 



^ = -2i (\u\ 2 - 2A 2 ) h - 2i {-u* x + 2i\u*) f (53) 

It is shown in Ref.(|l^]) that the condition on the initial function u {x, 0) for 
there to be only a finite number of nondegenerate points in the main spectrum 
is equivalent to the requirement that /, g and h be finite-order polynomials in 
the parameter A, which we take of degree N +1: f (x, t; A) = J^fJo 1 fa ( x > ^) ^ ) 
g (x, t; A) = £££ 9j (x, t) A J , h (x, t; A) = E^V h (x, t) X j . From the Eqs.(|| 
- 15^) it results however that g and h have degree N in A. 

One can check that the following combination is invariant in time and space 
and we note that it is actually the square of the Wronskian: 

|(/ 2 - 9 '>) = -|(/ 2 -»<■)=» (54) 

f -gh = -I [»'(«+, 

The Wronskian only depends on the spactral parameter A. We also know that for 
a subset of the main spectrum, the nondegenerate eigenvalues, the Wronskian is 
zero. Then the number of the nondegenerate eigenvalues is 2N + 2 (the degree 
of f 2 — gh as a polynomial in A) and we express the function f 2 — gh as a 
polinomial with constant coefficients (not depending on x and t). 



2JV+2 2N+2 

i [^(^ + ,^-)] 2 = / 2 - 5 ^ J P(A)- £P fc A fc = J] (A -A,) (55) 



fc=l j"=0 
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4.4.2 Product expansion of g and its zeros (ij (x,t). Introduction of 
the hyperelliptic functions ■ 

The functions g and h are both of order N in A. For g we will note the TV roots 
by ^ (x, t). The coefficient of X 2N+2 in f 2 - gh is f 2 N+1 . Due to Eq.(||) it is 
a constant which can be taken 1. Now since /jv+i = 1 the coefficient of X N in 
g is in (x,t). Written as a product 

N 

g = iu(x,t) JJ [X-fijix,^] (56) 
j=i 

By similar arguments we find that the coefficient of X N in h is m* (x, t) and 
the zeros of h are //* (a:, t). The functions /i^ are called the hyperelliptic 
functions. It will be proved below that finding the hyperelliptic functions leads 
immediately to the solution u (x, t). 

Now we calculate the expression in Eq. (|55|) for A = a zero of the function g, 
i.e. X is equal to hyperelliptic function jii (x,t): 

f 2 {x, t;X = /i m ) - gh = f 2 (x, t; fj, m ) = P (A = p, m ) 

or 



/ fa> *5 Mm) = a m V P (Mm) ( 57 ) 

Here the factor o~ m is a sheet index that indicates which sheet of the Riemann 
surface associated with yj P (A) the complex fx m lies on. 

Let us calculate from Eq.(f49|) and Eq.(|5^) the derivative at x of fi m : 

dg(x,t;X) ,du(x,t;X) -i-r , , 

N d[i (x t) N 
-iu(x,t;X)^2 — 3 -q^- II ( A_ MfcOM)) 

Replace here a zero of <?: A = /x m 

dg{x,t;fi m ) . dfi m (x,t) JX- 

— = -iu(x,t;fj, m ) 1- ]__[ (fi m - fj, k (x,t)) 

k—l^k^m 

On the other hand we have, from Eq.(^) 

^ = -2i^ m g(x,t;X = fi m ) -2u(x,t;X = n m ) f (x,t;X = fi m ) 
= -2uf 

then, since we have calculated / (x, i; A = Eq(p r 7|), we obtain 



dn m {x,t) -2io mX JP (fi m ) 



dx 



Tlk=l,kjtm (Mm — Mfc *)) 



(58) 
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for m = 1,2, ...,N. 

The coefficients of the term with A^ in the equation ( (49| ) are matched and 
we obtain iu x = — 2i<?at_i — 2it/jv .Using Eq.([56|) we find 

N 



d x In u = 2i [ E Mj + /. 



A' 



From the Eq.(p5|) we obtain the coefficient of the term A , i.e. /at 

2N+2 

2 



jv 



2N+2 



fc=l 



Then from the preceding equation it results 



(N ^ 2JV+2 \ 

j=i fe=i / 



(59) 



The same procedure is applied to the equation for c/ t , (|52|) and we obtain 
dn d (x,t) ( ^ 1 2 ^ 2 \ 9/^ (a:,*) 



St 



/ l 2 ^ 2 \ 

-2 E ^™ - 2 E A& 

\m#j fe=l / 



<9x 



(60) 



<9* In u — 2i 



-U 



3 / 2W+2 
E A J Afe _ 4 I E Afc 

j>k \ fc=l 

2N+2 \ / JV 

J E Afc ) E 



(61) 
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E 

j>k 



We can see that the problem has been reformulated: from the equation 
for the function u(x,t) the Lax-operators formulation leads to the problem for 
Bloch functions, ip + and ip~ and their Wronskian; then, using the "squared" 
eigenfunctions /, g and h we arrive at a formulation for the hyperelliptic func- 
tions fij (x, t) . Finding fij (x,t) leads immediately to u (x, t). 

These operations have a significative geometric counterpart: every point of 
the complex plane of the spectral parameter A is mapped, via Eq.([55"l) on the two- 
sheeted Riemann surface \J —\W (x, t; A) = yJP (A) with singularities at the 
nondegenerate points of the main spectrum, Xj. Removing the indeterminancy 
(by cutting and glueing the two sheets) we obtain a hyperelliptic Riemann 
surface of genus g — N. It results that we have to consider the variables 
fij (x,t), j — 1,N , as points on this surface but the motion of fij (x,t) is not 
simpler than the equation itself. 
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Now the next steps would be: (1) using the initial condition ^t.- (0, 0) on the 
function fj,j [x, t) we can solve the two equations ( p8| ) and (|6(]) and find fij (x, t); 
(2) then using the initial condition \u (0, 0)| for the function u (x, t) we can solve 
the equations (|5^) and (^Tj). The procedure will be: 

• Take the parameters j = 1, 2, 2N + 2} as known; these parameters 
are non-degenerate eigenvalues from the main spectrum, they can be 
found from the equation A (A) = ±2 and A is determined from u (x,0), 
the initial condition. 

• Choose initial conditions fij (0,0) and |u(0,0)| (see below) 

• solve the equations for (x,t) ; this means to find the hyperelliptic func- 
tions. Then solve the equations for u (x, t). 

4.5 Solution: the two-sheeted Riemann surface (hyperel- 
liptic genus- N Riemann surface) 

The variable \ij 's are points lying on the two-sheeted Riemann surface associated 
with 

1/2 



(2W+2 \ ' 

II (A-A fc )J =R(X) 



which has branch cuts at each of the nondegenerate points Afc. To 
eliminate the indeterminacy related to the square-root singularities (located 
at Aj's) the surfaces must be cut and reglued, obtaining a complex manifold 
of complex dimension one, a hyperelliptic Riemann surface. We need some 
constructions on this surface. 

4.5.1 Holomorphic differential forms and cycles on the Riemann sur- 
face 

On this hyperelliptic Riemann surface (denoted M) it is possible to define N 
linearly independent holomorphic (regular) differentials. The following is the 
canonical basis of differentials one-forms defined on the manifold M: 



R(X) 
XdX 



dU2 -R(xy 



All - Xl dX 

dUj - Rjxy 



:-S8 



A 



JV-l 



dX 



dU N EE 



R(X) 



On the surface M there are 27V topologically distinct closed curves (cycles). 
On a torus (TV = 1, called elliptic surface) there are two curves which cannot be 
deformed one into another. A more general Ricmann surface, with TV > 1 can 
be shown to be equivalent to a sphere with N handles and is called hyperelliptic 
Riemann surface. The number N is called the genus of the surface. We have 

N = genus of the surface = number of independent holomorphic differential forms 

There are 27V cycles which are split into two classes: aj cycles and bj cycles. 
Each of these cycles has a specified direction (an arrow) attached to it. To 
construct the cycles, the rules to be applied are: 

1. aj cycles do not cross any other aj cycle; bj cycles do not cross any other 



2. the cycle intersects bk only once and intersects no other b cycle; 

3. the intersection is such that, at the point of intersection, the vector tangent 
to the cycle au, the vector tangent to the cycle bk and the normal to the 
tangent plane of M represent a system of three vectors compatible with 
the orientation of M. 

This can be represented as: 



4.5.2 Periods and matrices of periods 

With the set of holomorphic differential forms and the set of cycles we can define 
the matrices of periods, A and B 



and the matrices A and B are invertible. A change of basis of holomorphic 
differential forms is a linear transformation represented by the matrix C: 



bj cycle; 



ooa = 0, bob = , djobk 






N 




fc=l 
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We can choose the matrix C as 



C = A^ 1 



Then the periods in the new basis becomes 




N N 




fe=i Ja " fe=i 



#j = Tj n 



b, 



where 



It can be shown that the matrix r is symmetric Tj-fc = T/y and has a positive- 
definite imaginary part Imr > 0. 

4.5.3 The Abel map 

The Abel map is defined from the Riemann surface M to the space C N and 
associates to a hyperelliptic Riemann surface (a sphere with TV hadles) a N- 
torus. Using the differentials dipj one constructs a change of variables by the 
following procedure: 

• choose a base point on M and call it p ; 

• define the variables Wj (x, t) as integrals on the Riemann A-surface of the 
differential forms from the arbitrary point po to the point which is the 
hyperelliptic function [i k (x,t): 



These variables have the remarkable property that their x and t depen- 
dence is trivial 





d 



N 



N 



to— 1 dn k 




and using 



dfj, k (x, t) 



= -2«(7 fe 



dx 



Un^k (Mfc - Mn) 
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it results 



N 



N 



dx 



m—1 



~"{ Un^k i^k ~ Mra) 



To calculate the sum 



JV 

E 



^r 1 



one can use the contour integral 

\ m ~ 1 d\ 



2™ lc n, 



^k (Pk - Mn) 



where the contour C encloses all of the poles n n counterclockwise. By the 
residue theorem 



N 

E 



Mr 1 



;s n„^ fc (MA - m„) 

One can evaluate this integral by compactifying the A plane into a sphere 
and noticing that the contour encloses the pole at z = oo. Evaluating the 
residuu at z = oo 



JV 

E 



vT 1 



= (5 



m,N 



and from this we obtain 



ffa W,-( a: ,t) = -2iC J -, JV = — kj 



A similar calculation gives 



d 
It 



-W,(x,t) 



-Ai 



c 



= — n, 

It results from these calculations 

W j (x,t) = ^-(k j x + Cl j t + d j ) 
where dj is a phase which is determined by the initial condition on fj, k . 

The Abel map has linearized the motion of the points fj, k (x,t). Now, in 
order to determine the function u (x, t) we must invert the Abel map, returning 
from the variables W to fj,. This is the Jacobi inversion problem. 
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4.5.4 The Jacobi inversion problem 

The return from the variables Wj (x, t) of the iV-torus back to the variables 
H k (x, t) of the Ricmann surface M (notice we will only need particular combi- 
nations of the functions /z) can be done in an exact way using the Riemann 9 
function. 

The argument of 9 is the N-tuple of complex numbers ~z G C N . 



oo 



9 (z\t) — ^ ••• ^ exp (2mm ■ z + Trim ■ r • to) 

mi- — oo mw^-co 

with the notations rn-z = X)tLi m kZk and m-T-m = X^fej=i T jk m j m k- Adding 
to the vector z a vector with a single nonzero element in the position k we obtain 
from the definition: 

9 (z + e k \r) = 9 (z\t) 

which means that the 9 -function has N real periods. Adding the full fc-th 
column of the matrix r, we obtain: 

6 (z + T k \r) = exp (-2iriz k - mr kk ) 9 (z\t) (62) 

where r kk is the diagonal element of the r matrix. 

There are two quantities which must be determined in order to have the 
explicit form of the solution. The space and respectively time equations for 

,2 



u(x,t) depend on the quantity J2j>k Vj^k = \ 
So we must determine the following combinations of the variables fi^ : 

N N 

J2^ m (x,t) and ^/j^Oct) 

m— 1 m— 1 

In order to calculate these two quantities, we shall start by introducing a 
series of functions , ipj (p) defined on the Riemann surface M: 

ipj (p) = [ dip j 

Jpa 

where p is a point on M , dtpj is the j -th holomorphic differential form 
defined on M and pa is an arbitrary fixed point on M. The functions ipj (p) ar e 
multiple valued since the contour is defined up to addition of any combination 
of the cycles on M. 

Now consider the function, F (p): 

F{p) = 9{iP (p)-K\r) 

where 9 is the iV-dimensional Ricmann theta function associated with the surface 
M; ip (p) is the TV-dimensional column of complex functions ipj (p); K is an 
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iV-dimcnsional column of complex numbers independent of the point p. The 
function F (p) is multivalued on the Ricmann surface M for the same reason as 
ip: moving the point p on the surface such as to turn around one of the b cycle 
and returning to the initial position, the functions t/jj will add elements of the r 
matrix. This will make the function F to change as imposed by the properties 
of the 9 function, shown in Eq. ( |62| ) . 

In order to render the function F single- valued, we replace its domain M by 
a new surface, obtained from M by dissecting it in a canonical fashion, along 
a basis of cycles. This new surface, M* is simply connected and has the bord 
composed of a number of arcs equal to AN. 

This operation is necessary in order to render the function F (p) entire and 
allow us define the inegral 



Io = ^-.[ dlnF(p) (63) 
2m JdM* 



around the contour of the surface M* . Applying Cauchy theorem, the integral is 
the number of zeroes of F (p) on the surface M*. This number is TV (Riemann). 

We now impose the condition that the N zeros of F (p) coincide with the 
N points (/x ■ (x, t) , Oj) . This fixes the values of the N complex numbers 
Kj {j = 1,...,N). By doing so the contour integral ( |63| ) calculated with the 
reziduum theorem will involve the variables fij . 

Calculation of the two combination of fi's We must remember that the 
complex A plane is covered by the two-sheeted Riemann surface whose compact 
version is M. Further this is mapped by Abel map onto the Jacobi iV-torus. To 
a variable on the Riemann surface M (formally also on M*), say p, corresponds 
a certain A (p). One defines the following integrals, which are proved to be real 
constants : 



h = ^-.j X(p)d\n[F{p)]=A 1 
27 " JdM* 



1 

2irl 



\ 2 {p)d\n[F{p)] = A 



2 



dM* 



They can be evaluated by the residuu theorem. By the choice of the constants 
K's, the zeroes of the of the function F (p) are located at fj,j. Then the residues 
are just the integrand (A) calculated in fi plus the residuu at the infinite, ±oo: 



N 

T i = Y] Mm + Re s [A (p) d\nF (p)] + Re s [A (p) d In F (p)] 

— ' A — >oo+ A — " — 



= 1 



N 



- Res [A 2 (p)dlnF(p)] + Re s_ [A 2 (p) dlnF (p)l 

A — >-oo+ A — >oo 
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The reason to write two residues at infinity is that A = oo is not a branching 
point which means that there are two points on the manifold M corresponding 
to A = oo, one on each sheet of the surface. We note the value of ip (p) when 
p is the point on M* corresponding to A — > oo ± . 
The result is 

N ■ a 

or 

N ■ a 

with 9 ± = F (r ± - if). The expression of 7 2 

^=&?-i^ln[F(r+-if)ir( r --iO]+i|ln 



F (r~ - g) 
F (r+ - K) 



F (r" - K) 
F (r+ - -FQ 



F (r- - if) 
F (r+ - if) 



Return to the equations for the function u With these expressions we 
come back to the two equations for the two partial derivatives of u. 



d , d F(r+-K) n . A 2 ^X\ 

— Inn = — In — )— -( + 2zAi - i > A„ 

9a; 9x F (r - K) — 



m— 1 



9, 9 , F (r+ - K) . 

— Intt = — In — ; + zconst 

dt dt F (r~ - K) 

The solution Let us note 

w = Q 

2JV+2 

k = 2A 1 - A ™ 

m— 1 

which are "external" frequency and wavelength. 
The solution is 

9(W-\t) 

u (x, t)=u (0, 0) exp (iko - iuot) g ^ w+ ^ 

where 

wf = ±(k jX + n j t + sf) 

The phases Sf are the part of r^ 1 — ifj which is independent of (x, t). 
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5 Stability of the envelope solutions 



In this Section we return to the equations obtained by multiple space and time 
scales analysis. Then the meaning of the variables <p, y, t is: first order (enve- 
lope) amplitude of the potential fluctuations in the ion turbulence and respec- 
tively large space and slow time variables. We focus on a single NSE in order to 
examine the properties of the stability of the solutions. The equation has the 
generic form 



and represents the equation for the envelope of a fast oscillation arising as 
solution of the barotropic equation. On infinite spatial domain this equation 
has soliton as solutions. On a periodic spatial domain (which is our case) the 
solutions can be cnoidal waves and they may be modulationally unstable. The 
instability can be examined in the neighbourhood of a solution by imposing 
a slight deviation and studying its time behaviour by classical perturbation 
expansion. This analysis and comparison with the geometric-algebraic results 
have been carried out in Ref . ( |l3|l ) . 

Consider (j> (y, t) is a solution and perturb it at the initial moment t = as 



The question is to find the relative behaviour of the two functions <fi (y, i) and 
<fi(y,t). It they depart exponentially for time close to the initial moment, then 
4> (y, t) is an unstable solution. Suppose the initial solution <f> (y, t) is the envelope 
of an exact plane wave i.e. it is constant in space. As a function of time it must 
have the form <f) (y, t) = pexp (2iqt) and we find 



To this solution we add a perturbation 

4> (y>t) = pexp (2iap t\ {1 + e \A\ exp (iky — iSlt) + A 2 exp (—iky + iQi)]} 

where A\ and A2 are real coefficients. Linearising the equation for small e we 
obtain: 

{(2p 3 cr +pd - peek 2 ) exp (iky - iilt) + 2p 3 aexp (-iky + iQt)\ Ax + 
{(2p 3 a-pQ* -pak 2 ) exp (-iky + idt) + 2p 3 crexp (iky - ittt)} A 2 
= 

Taking the complex conjugate we obtain another equation for the two real co- 
efficients A\ and A 2 . The condition of compatibility of this system of equations 




(64) 



(j)(y,t) =<t>(y,0) + eh(y) 



4>(y,t) =pexp(2iap 2 t) 



(65) 
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(the determinant equals zero) can only be fulfilled if the factor multiplying the 
function of y and t is allways zero. This gives us the dispersion relation 

CI = ±^k (ak 2 - 4p 2 a) 1/2 (66) 

which shows that for 

\k\<2\p\JZ 
V a 

(i.e. for long wavelengths) the two solutions (f>{y,t) and (f>(y,t) diverge expo- 
nentially in time and the function <f> (y, t) is modulationally unstable. From 
Eq.(p6[) we notice that when a (the Landau coefficient) < and a (the disper- 
sion coefficient) > the frequency Q is real and there is no linear instability. 
This remark is useful in connection with the numerical calculations described 
above, which connect the physical parameters with the coefficients (a, a) of the 
NSE. Sampling the space of the physical parameters we must retain those points 
where the linear instability is possible. In the above formula a , a and p are di- 
mensionless constants and the variables of the equation are normalised : time to 
to and length to L. The value of p is connected with the frequency of oscillation 
of the uniform poloidal envelope w^^* (s -1 ) by 



p = vl^S? (67) 

and the condition of linear instability is for the poloidal wavenumber m 



In the case of the following parameters: a = 0.399, a = 1.452, d = 2ira/L 

phys 
unif 



2tt (1/0.3) = 20. 944, a^fi = 10 5 (sec" 1 ) we have 



For the Nonlinear Schrodinger Equation it is possible to follow the perturbed 



solution much further in time in an exact manner. In Rcf.(|13|) it is shown that 
from a solution it is possible to return toward the initial time and to recover 
the results of the linear analysis. This is possible because the equation can be 
exactly solved for initial conditions starting arbitrary close to the envelope of 
the plane wave. The case of the plane wave (uniform envelope) is particular due 
to the simple form of the main spectrum, as will be discussed below. 

5.1 Modulation instability of the envelope of a plane wave 
5.1.1 General method 

The geometric-algebraic setting of Inverse Scattering Transform is the appro- 
priate framework for understanding the nature of the stability of the exact 
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solutions of the NSE. This is because the solution is strongly dependent on 
the topology of the hyperelliptic Ricmann curve which in turn is dependent on 
the initial function. A perturbation of the initial condition may change the 
main spectrum turning degenerate eigenvalues into new nondegenerate pairs of 
eigenvalues. This means that the hyperelliptic curve will have a different genus 
N' > N and new cycles and differential forms should be defined. The solution 
depends precisely on this topology. 

We order to investigate the stability of an exact TV-band solution of the NSE 
we have to develop in detail the construction of solution. In the particular case 
where the solution is the plane wave this analysis is simpler since all the spectral 
properties can be found analytically. It has been proved Jl|| that under the 
perturbation of the initial condition the degeneracies in the main spectrum will 
be removed and new degrees of freedom will become active. In the new solution 
some of the new degrees of freedom have little effect and can be neglected. 

5.1.2 Calculation of the main spectrum associated to the uniform 
amplitude solution 

The unperturbed solution is the envelope of the plane wave, Eq.(|65|). To de- 
teremine the main spectrum, the solution is considered initial condition (at 
t = 0) and introduced in the Lax operator. The eigenvalue problem of the Lax 
operator has the form 



where 




The equation is: 

a^ + ^( 1 + A )^ = o 

The general solution is written 

<f) 1 = A cos (ny) + B sin (/cy) 

v ° 

The second component will be also a combination of harmonic functions: 

4> 2 — C cos (ny) + D sin (ny) 
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The equation can have two independent one-column solutions. The one corre- 
sponding to the boundary conditions 



*i(» = 0;A) 



is determined by A = 1, B = ^Ja/aX/ (in) , C = and D = ^/a/ai/n. The 
solution is 

^;A)=( C0S ^7M^) ( 69 ) 

V Wan ' 

The second fundamental solution corresponds to the initial condition; 





$ 2 (y = 0;A)= ( J 



and it results: 



/a sm «y 



*»^ A )=(oos«y + Ai^=^ J (70) 

The fundamental solution matrix <f> is a 2 x 2 matrix having as column the 
two above solutions. From the fundamental solution matrix we calculate the 
monodromy matrix 

M(A) = $(y = d;A) 
The discriminant is the trace of the monodromy matrix. 

A (A) = 2 cos (d^Vl + X 2 ^j (71) 
We can proceed to the calculation of the spectrum. 



5.2 The spectrum of the plane wave solution 

5.2.1 The spectrum of the unperturbed initial condition (case of 
plane wave) 

We take a set of constants as derived from the numerical calculations that 
connects the physical conditions to the coefficients of the NSE, d = 2w (^) , a = 
0.399, a — 1.452. With these constants we calculate the spectrum on the 
complex A-plane. Although we can have analytical solutions the graphs below 
show the difficulty of resolving the details of the variation of the functions. To 
be noted that the function is intersected with a plane arbitrarly placed at the 
magnitude TrM = 5, for easier 3d-representation. 

This graph shows the presence of roots on the imaginary axis of A of the 
equation \TrM (A) — 2| =0 with very narrow variation of the function. Actually 
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Figure 14: Absolute value of Tr(M) — 2. The zeroes on the Im\ axis are difficult 
to resolve. 

it is easily to see that for the parameters adopted here, there are three roots of 
this equation, while in the graph we can hardly see two of them. A graphical 
better way of finding these roots is to look for the intersections of the contours 
of the two equations: 



where, again, the other two roots cannot be easily resolved. However the 
analytical expression allows to find the positions of the roots both on the Re A 
axis and on the Im A axis: 

By contrast, the variation along the real A axis easily evidences the zeroes 
of the corresponding to the main spectrum. 

The stability or instability of the Bloch fuctions (i.e. the values of the 
quantity m (A) which is the eigenvalue of the monodromy matrix) is governed 
by the discriminant A (A). We have —2 ^ A (A) ^5 2 for all real A and for 
A = ia with — 1 ^ a ^ 1. The main spectrum corresponds to the values of A for 
which A (A) = ±2. For real A, A = A" roal this means 



Re (TrM) -2 = 
Im [TrM) = 
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Along the \mX axis: lm(TrM)=0, Re(TrM) can have three zeros 




-0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

Mimag) 



Figure 16: 
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Figure 17: 



for n > 3. For pure imaginary A, A = iX° n 



A = + 

^n,imag 



1 - 



d\faJo 



1/2 



for rt = 0, 1,2 (with the values of the parameters a, a and d). These values are: 



±i, iA?, 



±i0.9422 and ^ 



±i0.7424. Since d is finite 



only a finite number of degenerate pairs will appear on the imaginary axis. At 
n = there are two eigenvalues Ao = ±z . These eigenvalues are nondegenerate 
in the sense that there is only one Bloch function for each of these values of A. 
This can be shown directly, constructing the Bloch functions form the matrix 
of fundamental solutions j69| ) and (70) with k — > (because A (A) = 2 implies 
K = 0) and showing that only one Bloch function can be obtained. 



5.2.2 The spectrum of the perturbed initial condition (case of plane 
wave) 

Starting from an initial condition slightly different of the plane wave 's envelope 



u(y,0) ->■ u(y,0)+eh (y) 



(72) 



we shall have to repeat all steps in the determination of the spectrum. This must 
be done in a comparative manner, identifying the degenerate eigenvalues Xj of 
the main spectrum which becomes nondegenerate after perturbation. The time 
behaviour of the exact solution developing from the perturbed initial condition 
is dominated by the imaginary part of the variables W which are obtained on 
the basis of the new (perturbed) quantities fj,j and this can be calculated. 
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|TrM-2| 




Figure 18: The graph of Tr(M) — 2 showing the zeroes on the ReX axis. 
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Using the new initial condition Eq.(|72|) we shall calculate the new main 
spectrum, i.e. the new values of the parameter A for which the discriminant 
is ±2. All the spectrum of the simple initial condition u (y, 0) is perturbed by 
order-e quantities. The main effect of this perturbation is that the degeneracies 
are broken. In order to construct the exact solution corresponding to the new 
position of the nondegenerate A values, we need to specify the cycles on the two- 
sheeted Riemann surface. The steps are as usual and consists of first choosing 
the cycles on the Riemann surface and the base of the holomorphic differentials; 
then, compute the matrices of periods (the matrix of a periods, and the matrix 
of the b periods); compute the t matrix; construct the Abel map, which maps 
the points of the compact two-sheeted Riemann surface onto the Jacobi torus; 
find the wavenumbers and the frequencies in the Jacobi manifold; make the 
Jacobi inversion, using the theta functions; 



Finding the base of the holomorphic differentials One of the eigen- 
value pair of the unperturbed system is Consider now that under 
the perturbation a degenerate eigenvalue A° splits into the conjugated pair 
(A" — Ek, A° + £fc). In Ref.(|ll|) it was constructed a new basis as a linear 
combination of the canonical basis: 

ajt 1 /1 , /^0\2 rim/j _ \n) d ^ , . at 

dU j = — Vl+(A,) m , for j =1,2,..., N 

where 

TV 



R(X) = Vl + \ 2 H [{X-X°k-e k ) (X-X° k +£ k )] 



1/2 



fc=l 



To understand the choice of basis it is useful to look at the limit e — > where 
the basis of holomorphic differential becomes 

KmdUj = -LVl + (A°) 2 - 1 -. dX nN 

This expression shows that, at the limit e — > 0, each differential "sees" only one 
pole, at Xj (which is a double point at this limit), and the square-root branch 
points ±i. 

The loop (cycle) aj encircles the pole A° and the period is 
lim Ak-i = lim / dUi- = / IvaxdUh = S 



ikj = nm / auk = I iimciujt = ojy 

This is because the cycle aj surrounds the pole A° or it does not surrounds any 
other pole. 

We find that, under the limit operation e — > 0, the matrix A goes over to 
the identity matrix to the first order in e, O (e) and that the matrix B becomes 

lim-Bjfc = T jk + O (e) 
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To see what happens with the B matrix in this limit, we must examine the 
off-diagonal terms and the diagonal terms. 

The off- diagonal terms are well-behaved and can be integrated analytically: 

Ym\Bjk = I \imdUk for i ' ^ k 
*-o Jb,^ 

and 

The contour bj passes through the middle of the segment connecting the two 
new Xj and goes around the pole at A = —i, for example. If this contour does 
not surround the pole at A° then the integral is elementary. If the contour 
intersecting the cut from -i to infinity encircles the pole at A° then it may 
be chnaged such that to intersect the other cut, connecting the branch point 
+i to infinity, and so it does not encircle any pole. In conclusion since the 
off-diagonal elements of the matrix r depend (to (e)) only of the positions of 
the original double points, which are determined only by the parameter d, it 
results that: the off-diagonal terms of the t matrix does not carry any 
information on the perturbed initial condition. Only the diagonal terms 
of the r matrix are affected by the initial condition. 

The diagonal terms of the r - matrix are singular at the limit e — > 0. 
Consider 

2m A,v/rT^( A -A° fe ) 2 -4 

Suppose that 

A° is on the imaginary axis 

Then we have 

2- V {J-i *JaJ VTT^^A-A") 2 -^ 



The integral can be made with elliptic functions but an approximation is (near 
£fc = 0) 

1 i 

T fefe = - - - In |e fe | + O (1) + O (e fc ) 

Z 7T 
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Consider now 



Af. is on real axis 

then the integral defining the diagonal elements of the matrix r (the integral 
along the b k cycles) is 




where the second integral is along the real axis. 

The first integral is done with the residuu theorem after taking safely the 
limit e k — and gives a finite imaginary contribution which is well-behaved for 
this limit. The second can be done and the result is 

T kk = --\n\e k \+iO(l) + 0(e k ) 

7T 

We then conclude that for e k — > the b - period matrix r becomes logarith- 
mically singular on the main diagonal. 

The main conclusion is: (1) The real-A axis degenaracies have J7j real so 
they are linearly stable. The imaginary axis degeneracies have Qj imaginary so 
they are linearly unstable. 

This conclusion can have important practical effects: in our case the de- 
generate eigenvalues on the imaginary A-axis split and induce instability of the 
plane wave solution, i.e. of the envelope. At the same time the solutio evolves 
toward one of the stable (soliton-like) solutions. Still much work is required to 
examine the stability of that solution to additional perturbation. For example, 
the main spectrum of the solution of the type sec h{y) can be obtained from 
the discriminant whose form is very complicated, with many purely imaginary 
(dangerous) eigenvalues. 

6 The growth of the perturbed solution as a 
source of poloidal asymmetry and spontaneous 
plasma spin-up 

6.1 The minimal condition for instability of the poloidal 
rotation 

We have found that the space-uniform solution (poloidal symmetry) is unstable 
and the turbulence will self-modulate according to one of the possible soliton- 
like solution. In a time picture, we can say that the poloidal symmetry will 
be left and the turbulence amplitude will develop a poloidal asymmetry which 
grows exponentially in time. The onset of the poloidal asymmetry will induce an 
asymmetric rate of transport on the poloidal direction p2[ . When the particle 
diffusion is poloidally asymmetric and the local particle confinement time is 
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Figure 19: The discriminant TrM (A) of the sec h (y) solution, to compare with 
the plane wave solution's discriminant. 



shorter than the damping time of the poloidal rotation, the poloidal rotation 
is unstable. The plasma can spontaneously spin-up. The growth rate of the 
asymmetry (which is given by the rate of departure from the uniform solution 
toward a soliton-like solution) determines the growth rate of the torque applied 
to plasma in the poloidal direction but the relation is not direct, depending 
essentially of the momentum balance in the poloidal direction. Any torque 
acting to generate a poloidal rotation must overcome the neoclassical magnetic 
pumping, which is an efficient damping mechanism. A balance of these forces 
depends on the comparison between the rates of (1) asymmetry growth, the 
radial variation of the asymmetry-induced terms, on one hand, and (2) rate of 
damping through magnetic pumping, on the other hand. 

For the simplest case of the space-uniform envelope solution, we have the 
instability has the linear growth rate 



,1/2 



where k < 2pyaja is here the wavenumber of the perturbation of the envelope. 
In the case that has served as example in the preceeding Section, we have 

p = to^Znif / (2°") ~ 1-86 and k < 7.08. Taking for example k — 3 we obtain 

7 e ™ = 7.67 or YpTys = 7 - 67 x 104 ( s_1 )- 



6.2 The soliton developing from the initial perturbation 

Since we know that the uniform state can be unstable we have to investigate 
the problem of the evolution of the system after a perturbation is applied and 
look for stable solutions. The solitons of this system can be realisable in real 
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conditions since they are intrinsically more stable. We will look for soliton-like 
solution for the uncoupled equations, taking a single NSE like (|64|). The time 
dependence can be 

<fi = exp (i/x 2 ai) w (x) 

giving the reduced equation 

r) 2 m „ (To 

fi z w + —w A = 



dx 2 ^ ' a 

where fi is a real constant. This equation has the following solution 

/2a 



-/isech (fix) 



The constant \x can be related to the modulation wavenumber whose upper 
bound is (^). Taking /i ~ 3 and a — 0.399 we obtain the frequency of steadily 
oscillating soliton u™ = /.i 2 a ~ 1.2, which is slower than the growth rate 
of the deformation, j lm — 7.67 obtained above (all in absolute units). Since 
the rate of change of the poloidal velocity (which will be calculated below) has 
comparable magnitude with w™ we conclude that the rotation has the time 
to develop before the oscillation reverses the sense of the applied torque. The 
maximum amplitude of the deformation is y/2a/afi = 2.224. Large smooth 
deformations (small /i) with less localised profile will have even slower steady 
oscillatory motion and however significant amplitude, and we can expect they 
are favored in real situations. 



6.3 The poloidal asymmetry of the diffusion fluxes 

Our approach, which is adapted to the investigation of the stability of the tur- 
bulence envelope, can only provide a model of the turbulent field, according to 
Eq. ( pl[ ) . This is composed of the primary modes propagating in the poloidal di- 
rection with amplitude modulated by the slowly varying envelopes obeying 
the system of coupled Nonlinear Schrodinger Equations ([50]) and (|3l|). In the 
radial direction the field has a variation given by the functions ip 1 2 (x) obtained 
by solving the equations (p2]). We will need to reconstitute the turbulent field 
in order to calculate the modified diffusion coefficient. We choose the following 
set of parameters: 

a — 1 (m) 
B T = 3.5 (T) 

T e = 1 (KeV) , T t = 1 (KeV) 



n = 10 20 (m~ 3 ) 
L n = 0.1 (m) 
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U = 6 x 10 5 (m/s) 



The extension on the radial direction of the mode is taken L = 0.1 (to) and 
this is also the normalising length. Other quantities are: normalising time 
to = 0.16 x 10" 6 (s) and n z = 0.335 x 10 9 (s^ 1 ), the ion-cyclotron fre- 
quency. From this it results the two parameters of the barotropic form of 
the ion mode, (3 = 111.75 and e = 0.476 x 10~ 2 . We take the following 
frequencies for the primary poloidal modes: lo\ — 0.1451 x 10 6 (s _1 ) and 
u>2 = 0.1931 x 10 6 (s _1 ). From the solution of the equations (|22| ) the two 
eigenvalues result: k\ — 43.51 (to -1 ) and k 2 = 56.81 (to -1 ). Solving the series 
of differential equations described in detail in Section 3 we obtain the following 
constants in the coupled Nonlinear Schrodinger Equations: a,\ = 0.1756 x 10~ 2 , 
a 2 = 0.2278 x 10~ 2 and respectively <j x = 0.18272 , a 2 = 0.5110. The frequency 
of the uniform oscillations, Eq.(|67|) is taken p = 0.0675 giving the maximum 
poloidal number for the large scale instability, Eq. (|68|) , to = [5.89] = 5. The 
growth rate of the modulational instability is 7 = 4175 (s _1 ). The expression 
of the field becomes 

4>{r,6) = <j) [(l+Ai(9)aech(y))cos(kiy-u)it + 8i)ip 1 (r) 
+ (1 + A 2 {0) sec h(y)) cos (k 2 y - uj 2 t + 5 2 ) ip 2 (r)] 

where <p is explained below. The initial phases 8\^ 2 are arbitrary and in the 
numerical calculations are taken as random quantities with uniform distribution 
on the interval (0, 1). 

We notice that the equations (p3) fulfilled by (p 1 2 are homogeneous and the 
amplitudes of the two functions are arbitrary. We fix the amplitude <p with 
the condition that the uniform profile of the averaged square of the potential 
of the fluctuating field obtains a typical value of the diffusion coefficient, Dq ~ 

|(/> | 2 ^}, which is taken D ~ 1 (to 2 /s). The following two figures represent the 
fluctuating potential, resulting from the combination of the factors described 
above. 

A simple way to take into account the effect of the turbulence amplitude 
modulation on the diffusion coefficient is to consider that locally the diffusion 

coefficient is proportional with the quantity ^|</>| 2 ^, where the averages should 
be taken over a statistical ensemble of realizations of the fluctuating field. This 
requires a carefull examination when applied to numerical calculations. We have 
to remember that the diffusion is a statistical process where a test particle inter- 
acting with fluctuating background field performs random displacements whose 
mean square grows linearly in time. The random displacements have to belong 
to a Gaussian statistical ensemble and this requires that all kind of elemen- 
tary experiences (i.e. displacements) to be represented with the corresponding 
probability. This is practically irrealisable in numerical calculations since the 
statistical ensemble of the stochastic field configurations is necessarly finite, i.e. 
incomplete as a Gaussian ensemble. Any average will be made on subensem- 
bles, more or less incomplete. The collection of field configurations are obtained 
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Figure 20: Instantaneous amplitude of the potential fluctuation, as calculated 
from the different contributions 
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Figure 21: Contour lines of the fluctuating potential; the radial extension is 
artificially extended to improve the picture. 
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radial 



Figure 22: The radial and poloidal dependence of the average diffusion coeffi- 
cient. 

from time and from space realizations of the field. The test particle should have 
the possibility to manifest the linear mean square displacement specific to the 
Gaussian ensembles and this rises the question on what is the minimal spatial 
extension on which the average should be done. We need in any case a space 
region larger than the correlation length of the fluctuating field which can be 
taken as the inverse of the characteristic wavenumber in the field (minimum of 
fci, k 2 ). 



The profile of the quantity D Q ) is represented in the figure below. The 

exemple chosed for these graphical representation is characterized by the im- 
posed constrain to obtain a rapid variation of the functions ip 1 and (p 2 on the 
radial r-interval, with the intention to simulate short wavelength ion turbulence, 
frequently observed in the numerical simulations of the Ion Temperature Gra- 
dient driven turbulence. This is achieved by asking the numerical procedures 
which solves the Schrodinger-type equations (^2|) to look for high order k of 
the eigenvalues. The pictures represented in this section correspond to K = 19 
for both eigenfunctions. This generates a potential dominated by small spatial 
scales and consequently, the fast variation of the local value of the diffusion 
coefficient. We take, in the case k = 19 a space region containing at least 
three complete oscillations of the typical wavelength. This should allow at least 
several random steps in the typical displacement of the particle. When we do 
this (taking arbitrarly a factor of few units larger than the correlation length) 
we note that the diffusion coefficient exhibits space variations. These are nat- 
ural and are inherited from the space dependence of the radial eigenfunctions 
ip 1 2 {r). Since in real situations the radial eigenfunctions are superposed with 
random initial phases we need to enlarge the region of averaging over most of the 
radial interval considered, up to a level of saturation of these spatial variation 
which can be considered convenient for the estimation of the poloidal torque. 
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Figure 23: Radial variation of the diffusion coefficient after poloidal averaging. 



6.4 The poloidal torque applied on plasma 

We use the equations describing the plasma poloidal and toroidal rotation de- 
rived in 0, U): 

dn 1 d , _ . 
— + - — (rnvr) = 
at r or 



Id 

\nV v ) + - g^ [rn (y v v r - qVev r )] = 



d dV d 

[V v + 6 (l + 2q 2 ) Vg] + Vr-Q-^- ~ v r-Q— {oY v ) + magnetic pumping = 



The definition of the variables is n = (n (r)), nV v (r) = (^nv v ^j and Vg (r) — 

vg-^^i where v v and vg are local toroidal and poloidal plasma velocities in 
the magnetic surface. The magnetic field is 

with R = Rq + r cos 9 and Q= j^^- surface average is made according to 
the formula 
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In the collisional regime the magnetic pumping term determined by parallel 
viscosity is less effective and the spin up can be expected. The decay rate of 
the poloidal rotation by this mechanism is pTS(] : 



3 / 1 \ 1 / I N " 



^I^+vJ (?3) 

where / is the mean-free path. The velocities v r and v r arise from the diffusive 
flux in the radial direction, to which is associated the local radial velocity V r . 
We define the two velocities 

V r = (V r ) 

v r = (2cos6> V r ) (74) 



If there is no poloidal variation of the local radial velocity V r the average in ( |74D 
is that of the trigonometric function cos# and the result is zero. The average 
will not be zero if there is a poloidal variation of the diffusion-generated particle 
radial velocity, V r (9). We can introduce the modified diffusion coefficient D 
related to the nonsymmetry, 

— — On 

nv r = r r = -D— (75) 
or 

with 

D = (d (\(/}\ 2 ) 2cos6>\ (76) 

\ \ / st I pol 

where the first averaging is statistical and the second is on the poloidal surface. 
The equation for the plasma poloidal velocity is 

© (1 + 2? 2 ) + iMpVe) + qVe^ (nrv r ) = (77) 

The average velocity Eq.(|7|) still dependes on the radial coordinate via D a , the 
density n (r) and its derivatives and, more important, via the space dependence 
of the potential cp (x) left after the statistical averaging needed for the calculation 
of the diffusion coefficient (|7"rj|). The following quantity must be numerically 
calculated 



and this directly provides the rate of change of the poloidal velocity due to the 
asymmetry. This must be compared to the magnetic pumping term ([73]) and 
we expect rotation instability when 

7 ~ Imp > 
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Figure 24: Radial variation of the rate of time variation 7 applied on plasma in 
the poloidal direction, due to the asymmetry. 

In this formula the soliton profile of the poloidal asymmetry is considered satu- 
rated and the motion of the solitons (which would give periodic variation of the 
torque) is neglected. As is shown in Rcf. |]l7f the high positive second derivative 
of the density profile, combined with a dependence of the diffusion coefficient 
Dq on the shear of the poloidal velocity are the favorable circumstances for 
the asymmetry-induced torque to overcome the magnetic pumping term and a 
spontaneous spin up to occur. We do not take into account the variation of Dq 
with the velocity shear, but notice that the space variation of ( |76| ) also con- 
tribute to obtain a positive value for 7 sufficient for being comparable to Jmp 
and even to be larger. The most important contribution comes however from 
the positive- valued second density gradient, which can be seen as an indication 
that the spin-up can be expected in such region on the minor radius. 

The values represented in the figure have been calculated using a gaussian 
profile for the density. The decay due to magnetic pumping is of the order of 
Imp ~ 12 (s _1 ). We note that the torque applied to plasma shows a radial 
variation essentially inherited from the space variation of the statistically aver- 
aged potential fluctuations, which in general are not monotone functions of the 
radius on time scales of several eddy turn-over times. These spatial variations 
of the diffusion coefficient generate the profile of the torque which change the 
shear of the seed velocity and provide a characteristic which is expected for 
zonal flows. Even outside the domain of the parameters where the spin up is 
efficient compared to the magnetic pumping, the soliton asymmetry and the 
Stringer torque can contribute to the action of other mechanisms to generate 
poloidal rotation. 



63 



7 Discussion 



We have examined the problem of stability of the envelope of the ion wave turbu- 
lence. Starting from a simple model for the ion dynamics and retaining the full 
nonlinear contribution of the ion polarization current we have derived an equa- 
tion of the same form as the barotropic model of the atmosphere. The multiple 
space-time scale analysis shows that the envelope of the turbulence results from 
a set of coupled cubic nonlinear Schrodingcr equations. The simplest solution is 
the space-constant envelope corresponding to a symmetric average amplitude of 
the turbulence on the poloidal circumference. In order to examine the stability 
of this solution we first review the method of inverse scattering transform in 
the geometric-algebraic framework, able to obtain the exact solution developing 
from a certain class of initial conditions (finite-band potentials) . In this frame- 
work appears clear the origin of the sensitivity to the perturbations, connected 
to the topology of a hypcrclliptic Ricmann surface. The uniform solution is 
unstable and evolves to one of the soliton-type solution. The nonuniformity of 
the turbulence average profile generates a torque on the plasma, which may be 
sufficient to overcome the poloidal damping. 

The simple physical model we have employed may be partly justified by the 
observation that the analysis is essentially independent of the other possible 
sources of poloidal asymmetry: the neoclassical variation of parameters on the 
magnetic surface, the toroidal effect on the nonlincarity represented as mode 
coupling, the effect of the particle drifts on the linear and renormalized propa- 
gators, etc. The cubic nonlinear Schrodingcr equation has many solutions on the 
periodic domain and the relevance for the self-modulation of the envelope of the 
turbulence should be examined. Here we only considered the uniform envelope 
since its instability can be sufficient to induce plasma poloidal rotation. This 
rotation has an origin different of the usual sources of rotation: neoclassical, 
ion-loss, Reynolds stress or momentum injection. Our approach is essentially 
mathematical and so it operates at a more formal level. The physical origin 
of the instability seems to be the tendency of the system (averaged amplitude 
of the turbulence) to evolve to profiles characterised by the balance of the dis- 
persion and the nonlincarity. What remains to be examined is the quantitative 
constraint on the minimal initial perturbation needed to produce a particular 
evolution: for a weak perturbation the solution may be simple "radiation" while 
for initially larger amplitude a soliton emerges ensuring higher robustness and 
a localised form. It is in this case that the asymmetry can induce plasma ro- 
tation. This question can be answered by examining the exact solution, with 
specific numerical tools for the determination of the main spectrum. It might 
be possible that the Riemann curve is not a hypcrclliptic surface but a manifold 
close to a finite-genus curve. 
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